function var=Equi(t,y,rho,g,n,frisch,kL,theta,beta,phi,pix,kW,RefWealth,sigma,NaturalRate,piW)
var=[y(1)*(max(NaturalRate+pix+phi*(y(2)-pix+beta*(kL*y(1)^(1+1/frisch)-1)),0)-max(y(2)+beta*(kL*y(1)^(1+1/frisch)-1),piW-g)-g-rho+kW*(y(3)-RefWealth)^(-sigma)*y(1)) %Euler equation
     theta*(max(y(2)+beta*(kL*y(1)^(1+1/frisch)-1),piW-g)-y(2)) %Inflation anchor dynamics
     (max(NaturalRate+pix+phi*(y(2)-pix+beta*(kL*y(1)^(1+1/frisch)-1)),0)-max(y(2)+beta*(kL*y(1)^(1+1/frisch)-1),piW-g)-n-g)*y(3) %Ponzi dynamics
     max(NaturalRate+pix+phi*(y(2)-pix+beta*(kL*y(1)^(1+1/frisch)-1)),0)]; %Nominal interest rate (y(4) is only used to compute the fall in bond prices induced by the path of the nominal interest rate)